## ---- Data ----
## SOEP data
load("dat/proc-data/soep_imp.RData")

## ---- Shortened TSCS: AfD analyses, 2014-2018 (demeaning) ----
## Variable selection for demeaning
current_dm_vars <- soep_imp$imputations$imp1 %>%
  dplyr::select(dplyr::ends_with("_cwu")) %>%
  names() %>%
  substr(., 1, nchar(.) - 4)
additional_dm_vars <- c(
  "nbh_met_p_deutschl",
  "nbh_alq_p_quote",
  "kr_uemprate",
  "kr_foreigner"
)
dm_vars <- c(current_dm_vars, additional_dm_vars)

## Updated regional data
regional_data <- rio::import(
  "../../soep-data/soep.v38/stata/regionl.dta"
) %>%
  dplyr::filter(kkz > 0) %>%
  dplyr::select(kkz,
                syear,
                kr_uemprate,
                kr_foreigner) %>%
  dplyr::rename(year = syear) %>%
  dplyr::mutate_at(
    .vars = vars(kr_uemprate,
                 kr_foreigner),
    .funs = ~ifelse(. < 0, NA_real_, .)
  ) %>%
  dplyr::distinct() %>%
  dplyr::group_by(kkz, year) %>% 
  dplyr::summarize_all(mean) %>%
  dplyr::ungroup()

## Data management
for (m in seq_along(soep_imp$imputations)) {
  soep_imp$imputations[[m]] <- soep_imp$imputations[[m]] %>%
    dplyr::mutate(yrs_at_current_address = year - hh_at_current_address,
                  age_at_current_address = age - yrs_at_current_address) %>%
    dplyr::mutate(ltr_since_3 = as.numeric(yrs_at_current_address >= 3),
                  ltr_since_5 = as.numeric(yrs_at_current_address >= 5),
                  ltr_since_7 = as.numeric(yrs_at_current_address >= 7),
                  ltr_age_30 = as.numeric(age_at_current_address <= 30)) %>%
    # filter(social_housing == 0 | is.na(social_housing)) %>%
    dplyr::mutate(lw = as.numeric(lef == 1 | spd == 1 | gre == 1),
                  rw = as.numeric(cdu == 1 | fdp == 1 | afd == 1)) %>%
    filter(renter_no_rent == 0) %>%
    dplyr::select(
      id,
      hh_id,
      plz,
      kkz,
      year,
      year_fac,
      weight,
      starts_with("cold_rent_load"),
      starts_with("cold_rent_sqm"),
      wr_ecego,
      wr_immig,
      afd,
      vote_afd,
      gre,
      lef,
      fdp,
      cdu,
      spd,
      non,
      lw,
      rw,
      owner,
      all_of(dm_vars),
      starts_with("rent"),
      starts_with("cmr_arm"),
      starts_with("log_hinc_eq"),
      starts_with("gtyp"),
      starts_with("prop_personal_hinc"),
      starts_with("hh_prop_ecact"),
      starts_with("hh_mmb"),
      starts_with("mover"),
      starts_with("lm_part"),
      starts_with("edu"),
      starts_with("age"),
      starts_with("east"),
      starts_with("fem"),
      starts_with("east"),
      starts_with("ltr_"),
      social_housing
    ) %>%
    dplyr::filter(not(is.na(rent_umn))) %>%
    dplyr::filter(not(is.na(rent_cwu))) %>%
    dplyr::filter(not(is.na(cmr_arm_umn))) %>%
    dplyr::filter(not(is.na(cmr_arm_cwu)))
}

## Local market data
load("dat/proc-data/plz5_f_und_b.RData")
plz5_f_und_b <- plz5_f_und_b %>%
  st_set_geometry(NULL) %>%
  filter(year %in% 2005:2018) %>%
  dplyr::select(plz, year, cmr_arm, rent) %>%
  distinct() %>%
  group_by(plz) %>%
  mutate(
    cmr_arm_2005 = mean(cmr_arm[year %in% 2005:2007], na.rm = TRUE),
    cmr_arm_2018 = mean(cmr_arm[year %in% 2016:2018], na.rm = TRUE),
    cmr_arm_chg = cmr_arm_2018 - cmr_arm_2005
  ) %>%
  mutate(
    rent_2005 = mean(rent[year %in% 2005:2007], na.rm = TRUE),
    rent_2018 = mean(rent[year %in% 2016:2018], na.rm = TRUE),
    rent_chg = rent_2018 - rent_2005,
    cmr_arm_chg = ifelse(is.na(cmr_arm_chg), rent_chg, cmr_arm_chg),
    cmr_arm = ifelse(is.na(cmr_arm), rent, cmr_arm)
  ) %>%
  ungroup() %>%
  group_by(year) %>%
  mutate(
    rent_cat =
      ifelse(
        year == 2004,
        NA_integer_,
        cut(
          rent,
          breaks = quantile(rent, c(0, 1 / 3, 2 / 3, 1), na.rm = TRUE),
          labels = paste0("T", 1:3)
        )
      ),
    cmr_arm_cat =
      cut(
        cmr_arm,
        breaks = quantile(cmr_arm, c(0, 1 / 3, 2 / 3, 1), na.rm = TRUE),
        labels = paste0("T", 1:3)
      )
  ) %>% 
  ungroup() %>%
  dplyr::select(plz, cmr_arm_chg, rent_chg, rent_cat, cmr_arm_cat, year) %>%
  distinct() %>%
  na.omit() %>%
  mutate(
    cmr_arm_chg_cat =
      cut(
        cmr_arm_chg,
        breaks = quantile(cmr_arm_chg, c(0, 1 / 3, 2 / 3, 1)),
        labels = paste0("T", 1:3)
      ),
    rent_chg_cat =
      cut(
        rent_chg,
        breaks = quantile(rent_chg, c(0, 1 / 3, 2 / 3, 1)),
        labels = paste0("T", 1:3)
      )
  )

## Additional data management
for (m in seq_along(soep_imp$imputations)) {
  soep_imp$imputations[[m]] <- soep_imp$imputations[[m]] %>%
    # top code household rents per sqm
    mutate(
      cold_rent_sqm = case_when(
        cold_rent_sqm > 40 ~ 40,
        TRUE ~ cold_rent_sqm
      )
    ) %>% 
    mutate(
      gtyp4 = case_when(
        gtyp %in% c(
          "[1] Kernstaedte>==gr.Verdraum",
          "[2] Kernstaedte<==gr.Verdraum",
          "[9] Kernstaedte<==Verdansatz"
        ) ~ "urban",
        gtyp %in% c(
          "[7] Mi-zent.=laendl=gr.Verdraum",
          "[10] Mi-zent.=verd.=Verdansatz",
          "[12] Mi-zent.=laendl=Verdansatz",
          "[14] Mi-zent.=verd.=laendl.",
          "[16] Mi-zent.=laendl=laendl."
        ) ~ "towns",
        gtyp %in% c(
          "[3] Mi-zent.=hochverd.=gr.Verdraum",
          "[5] Mi-zent.=verd.=gr.Verdraum",
          "[4] so.Gem.=hochverd.=gr.Verdraum",
          "[6] so.Gem.=verd.=gr.Verdraum"
        ) ~ "suburban",
        gtyp %in% c(
          "[8] so.Gem.=laendl=gr.Verdraum",
          "[11] so.Gem.=verd.=Verdansatz",
          "[13] so.Gem.=laendl=Verdansatz",
          "[15] so.Gem.=verd.=laendl.",
          "[17] so.Gem.=laendl=laendl."
        ) ~ "rural",
        TRUE ~ NA_character_
      ),
      gtyp3 = ifelse(gtyp4 == "towns", "rural", gtyp4)
    ) %>%
    left_join(plz5_f_und_b,
              by = c("plz", "year")) %>%
    dplyr::filter(weight > 0)
  
  ## Repair names
  names(soep_imp$imputations[[m]]) <- 
    gsub(" ", "_", names(soep_imp$imputations[[m]]))
}

## Join county-level data, within-demeaning
for (m in seq_along(soep_imp$imputations)) {
  soep_imp$imputations[[m]] <- soep_imp$imputations[[m]] %>%
    dplyr::select(-dplyr::starts_with("kr_")) %>%
    dplyr::left_join(regional_data,
                     by = c("kkz", "year")) %>%
    dplyr::select(-ends_with("_cwu"), -ends_with("_cwu")) 
}

# 2014:2018 Data for AfD support analyses
# Make time-invariant categorical moderators
soep_afd <- soep_imp
soep_afd$imputations <-
  lapply(soep_afd$imputations, function (imp) {
    imp %>%
      dplyr::filter(year >= 2014) %>%
      dplyr::select(-dplyr::ends_with("_cwu")) %>%
      dplyr::select(-dplyr::ends_with("_umn")) %>%
      group_by(id) %>%
      mutate(log_hinc_eq_umn = mean(log_hinc_eq, na.rm = T)) %>%
      mutate(cmr_arm_cat = get_mode(cmr_arm_cat)[1],
             cmr_arm_chg_cat = get_mode(cmr_arm_chg_cat)[1]) %>%
      ungroup() %>%
      mutate(
        log_hinc_eq_cat =
          cut(
            log_hinc_eq_umn,
            breaks = quantile(log_hinc_eq_umn, c(0, 1 / 3, 2 / 3, 1)),
            labels = paste0("T", 1:3)
          )
      ) %>%
      dplyr::group_by(year) %>%
      dplyr::mutate(weight = sqrt(weight / mean(weight))) %>%
      ungroup()
  })

## [2014, 2018] Data for vote choice analyses
soep_vot <- soep_imp
soep_vot$imputations <-
  lapply(soep_afd$imputations, function (imp) {
    imp %>%
      dplyr::filter(year %in% c(2014, 2018)) %>%
      dplyr::select(-dplyr::ends_with("_cwu")) %>%
      dplyr::select(-dplyr::ends_with("_umn")) %>%
      group_by(id) %>%
      mutate(log_hinc_eq_umn = mean(log_hinc_eq, na.rm = T)) %>%
      ungroup() %>%
      dplyr::group_by(year) %>%
      dplyr::mutate(weight = sqrt(weight / mean(weight))) %>%
      ungroup()
  })

## ---- Models ----
## Outcomes
outcomes <- c(
  "afd",
  "vote_afd"
)

## Main predictors
predictors <- predictors_aux <- c("cmr_arm")

## Renter controls
rent_vars <- rent_vars_aux <- c(NA_character_,
                                "cold_rent_sqm")

## Moderators
individual_moderators <-
  individual_moderators_aux <- c(NA_character_,
                                 "log_hinc_eq_umn",
                                 "log_hinc_eq_cat")
regional_moderators <- c(
  NA_character_,
  "gtyp3",
  "[_pred_]_cat",
  "[_pred_]_chg_cat"
)
regional_moderators_aux <- c(
  NA_character_,
  "gtyp3"
)

## Subsets
subset_owner <- subset_owner_aux <- c("owner == 0",
                                      "owner == 1")
subset_ltr <- c(NA_character_,
                "ltr_since_3 == 0",
                "ltr_since_5 == 0",
                "ltr_since_3 == 1",
                "ltr_since_5 == 1")
subset_ltr_aux <- "ltr_since_5 == 1"

## Model formula components
predictor_terms <- predictor_terms_aux <- 
  c("[_pred_]")

renter_control_terms <- renter_control_terms_aux <- 
  c("[_rent_var_]")

common_control_terms <-
  list(
    "none_year_fe" = "year_fac",
    "full" = c(
      "log_hinc_eq",
      "prop_personal_hinc",
      "hh_prop_ecact",
      "hh_mmb",
      "lm_part_Atypical",
      "lm_part_Inactive",
      "lm_part_Unemployed",
      "lm_part_In_Education",
      "lm_part_Pensioners",
      "edu5",
      "year_fac"
    ),
    "full_plus2" = c(
      "log_hinc_eq",
      "prop_personal_hinc",
      "hh_prop_ecact",
      "hh_mmb",
      "lm_part_Atypical",
      "lm_part_Inactive",
      "lm_part_Unemployed",
      "lm_part_In_Education",
      "lm_part_Pensioners",
      "edu5",
      "nbh_met_p_deutschl",
      "nbh_alq_p_quote",
      "year_fac"
    )
  )

common_control_terms_aux <- common_control_terms["full"]

## ---- Compile model code ----
model_formula <-
  'mod <- fixest::feols(
  [_yvar_] ~
  [_xvars_] | id,
  cluster = ~ id + plz,
  data = dat[[m]],
  combine.quick = FALSE
)'

model_formula_list_main <- list()
counter <- 0L
for (outcome in outcomes) {
  if (outcome %in% c("cold_rent_load",
                     "cold_rent_sqm")) {
    subset_owner_local <- subset_owner[1]
  } else {
    subset_owner_local <- subset_owner
  }
  for (predictor in predictors) {
    for (moderator in individual_moderators) {
      for (regional in regional_moderators) {
        for (owner in subset_owner_local) {
          if (owner == "owner == 1" |
              outcome %in% c("cold_rent_load",
                             "cold_rent_sqm")) {
            rent_vars_local <- NA_character_
          } else {
            rent_vars_local <- rent_vars
          }
          
          for (ltr in subset_ltr) {
            for (rent_var in rent_vars_local) {
              for (controls in names(common_control_terms)) {
                ## Counter
                counter <- counter + 1L
                
                ## Covariates
                if (is.na(moderator) & is.na(regional)) {
                  xvars <-
                    collapse_plus(predictor_terms) %>%
                    gsub("\\[\\_pred\\_\\]",
                         predictor,
                         .)
                  
                  if (not(controls == "none")) {
                    xvars <- paste0(xvars,
                                    collapse_plus(common_control_terms[[controls]]))
                  }
                  
                  if (not(is.na(rent_var))) {
                    renter_controls_local <-
                      gsub(
                        "\\[\\_rent\\_var\\_\\]",
                        rent_var,
                        collapse_plus(renter_control_terms)
                      )
                    xvars <- paste0(xvars, renter_controls_local)
                  }
                } else if (not(is.na(moderator)) & is.na(regional)) {
                  xvars <-
                    sapply(predictor_terms, collapse_by, moderator) %>%
                    collapse_plus() %>%
                    gsub("\\[\\_pred\\_\\]",
                         predictor,
                         .)
                  
                  if (not(controls == "none")) {
                    xvars <- paste0(xvars,
                                    collapse_plus(common_control_terms[[controls]]))
                  }
                  
                  if (not(is.na(rent_var))) {
                    renter_controls_local <-
                      sapply(renter_control_terms,
                             collapse_by,
                             moderator) %>%
                      collapse_plus() %>%
                      gsub("\\[\\_rent\\_var\\_\\]",
                           rent_var,
                           .)
                    xvars <- paste0(xvars, renter_controls_local)
                  }
                } else if (is.na(moderator) & not(is.na(regional))) {
                  xvars <-
                    sapply(predictor_terms, collapse_by, regional) %>%
                    collapse_plus() %>%
                    gsub("\\[\\_pred\\_\\]",
                         predictor,
                         .)
                  
                  if (not(controls == "none")) {
                    xvars <- paste0(xvars,
                                    collapse_plus(common_control_terms[[controls]]))
                  }
                  
                  if (not(is.na(rent_var))) {
                    renter_controls_local <-
                      sapply(renter_control_terms,
                             collapse_by,
                             regional) %>%
                      collapse_plus() %>%
                      gsub("\\[\\_rent\\_var\\_\\]",
                           rent_var,
                           .) %>%
                      gsub("\\[\\_pred\\_\\]",
                           predictor,
                           .)
                    xvars <- paste0(xvars, renter_controls_local)
                  }
                } else if (not(is.na(moderator)) &
                           not(is.na(regional))) {
                  xvars <-
                    sapply(predictor_terms,
                           collapse_by,
                           collapse_by(moderator, regional)) %>%
                    collapse_plus() %>%
                    gsub("\\[\\_pred\\_\\]",
                         predictor,
                         .)
                  
                  if (not(controls == "none")) {
                    xvars <- paste0(xvars,
                                    collapse_plus(common_control_terms[[controls]]))
                  }
                  
                  if (not(is.na(rent_var))) {
                    renter_controls_local <-
                      sapply(
                        renter_control_terms,
                        collapse_by,
                        collapse_by(moderator, regional)
                      ) %>%
                      collapse_plus() %>%
                      gsub("\\[\\_rent\\_var\\_\\]",
                           rent_var,
                           .) %>%
                      gsub("\\[\\_pred\\_\\]",
                           predictor,
                           .)
                    
                    xvars <- paste0(xvars, renter_controls_local)
                  }
                }
                
                xvars <- substr(xvars, 1, nchar(xvars) - 7)
                
                cat(xvars)
                cat("\n")
                
                ## Data
                if (outcome == "afd") {
                  data <- "soep_afd$imputations[[m]]"
                  
                } else if (outcome == "vote_afd") {
                  data <- "soep_vot$imputations[[m]]"
                } else {
                  data <- "soep_imp$imputations[[m]]"
                }

                if (not(is.na(owner))) {
                  data <- paste(data,
                                paste0("filter(", owner, ")"),
                                sep = " %>% \n    ")
                }
                
                if (not(is.na(ltr))) {
                  data <- paste(data,
                                paste0("filter(", ltr, ")"),
                                sep = " %>% \n    ")
                }

                ## Combine
                model_formula_local <- model_formula %>%
                  gsub("\\[\\_yvar\\_\\]", outcome, .) %>%
                  gsub("\\[\\_xvars\\_\\]",
                       xvars,
                       .) %>%
                  gsub("\\[\\_data\\_\\]",
                       data,
                       .)

                cat(model_formula_local)
                cat("\n")
                
                ## Return
                model_formula_list_main[[counter]] <-
                  list(
                    outcome = outcome,
                    predictor = predictor,
                    moderator = moderator,
                    regional = gsub("\\[\\_pred\\_\\]",
                                    predictor,
                                    regional),
                    owner = owner,
                    subset = ltr,
                    rent_var = rent_var,
                    controls = controls,
                    data = data,
                    formula = model_formula_local
                  )
              }
            }
          }
        }
      }
    }
  }
}

## ---- Estimation (main analyses) ----
dir.create("est/stayer_fe_analyses")

focal_analyses <-
  c(
    209,
    224
  )

## Capture start time
analysis <- 1
ptm <- proc.time()[3]
#for (analysis in analysis:length(model_formula_list_main)) {
for (analysis in focal_analyses) {
  ## Extract information for prediction
  mfx_vars <- 
    na.omit(unlist(model_formula_list_main[[analysis]][c("predictor", "rent_var")]))
  mod_var <- na.omit(model_formula_list_main[[analysis]]$moderator)
  reg_var <- na.omit(model_formula_list_main[[analysis]]$regional)
  
  ## Estimate model in parallel
  current_formula <- model_formula_list_main[[analysis]]$formula
  
  ## Data
  dat <- lapply(seq_along(soep_imp$imputations),
                function (m) {
                  eval(parse(text = model_formula_list_main[[analysis]]$data)) %>%
                    droplevels()
                })
  
  ## Model estimation
  est <- lapply(seq_along(soep_imp$imputations),
                function (m) {
                  eval(parse(text = current_formula))
                })
  
  ## Descriptives
  get_descriptives(
    dat[[1]],
    mfx_vars = mfx_vars,
    mod_var = mod_var,
    reg_var = reg_var,
    num_breaks = 31L
  ) %>%
    write.csv(paste0(
      "est/stayer_fe_analyses/",
      "descriptives",
      analysis,
      "_main.csv"
    ))
  
  ## General and variables-specific prediction values
  xmeans <- bind_cols(
    dat[[1]] %>%
      summarize_if(is.numeric, median, na.rm = TRUE),
    dat[[1]] %>%
      mutate_if(is.character, as.factor) %>%
      summarize_if(is.factor, get_mode)
  )
  if ("log_hinc_eq_umn" %in% names(dat[[1]]))
    log_hinc_eq_umn <- seq(
      quantile(dat[[1]]$log_hinc_eq_umn, .0005, na.rm = TRUE),
      quantile(dat[[1]]$log_hinc_eq_umn, .9995, na.rm = TRUE),
      length.out = 21L
    )
  if ("log_hinc_eq_cat" %in% names(dat[[1]]))
    log_hinc_eq_cat <- levels(dat[[1]]$log_hinc_eq_cat)
  if ("gtyp3" %in% names(dat[[1]]))
    gtyp3 <- levels(as.factor(dat[[1]]$gtyp3))
  if ("rent_cat" %in% names(dat[[1]]))
    rent_cat <- levels(dat[[1]]$rent_cat)
  if ("rent_chg_cat" %in% names(dat[[1]]))
    rent_chg_cat <- levels(dat[[1]]$rent_chg_cat)
  if ("cmr_arm_cat" %in% names(dat[[1]]))
    cmr_arm_cat <- levels(dat[[1]]$cmr_arm_cat)
  if ("cmr_arm_chg_cat" %in% names(dat[[1]]))
    cmr_arm_chg_cat <- levels(dat[[1]]$cmr_arm_chg_cat)
  
  ## Prediction data
  if (length(mod_var) == 0 & length(reg_var) == 0) {
    at_values <- NULL
    pred_data <- NULL
  } else if (length(mod_var) > 0 & length(reg_var) == 0) {
    at_values <- list(get(eval(mod_var)))
    names(at_values) <- mod_var
    pred_data <-
      eval(parse(
        text = paste0(
          "marginaleffects::datagrid(model = est[[m]], ",
          mod_var,
          " = at_values[[mod_var]])"
        )
      ))
  } else if (length(mod_var) == 0 & length(reg_var) > 0) {
    at_values <- list(get(eval(reg_var)))
    names(at_values) <- reg_var
    pred_data <-
      eval(parse(
        text = paste0(
          "marginaleffects::datagrid(model = est[[m]], ",
          reg_var,
          " = at_values[[reg_var]])"
        )
      ))
  } else if (length(mod_var) > 0 & length(reg_var) > 0) {
    at_values <- list(get(eval(mod_var)),
                      get(eval(reg_var)))
    names(at_values) <- c(mod_var, reg_var)
    pred_data <-
      eval(parse(
        text = paste0(
          "marginaleffects::datagrid(model = est[[m]], ",
          mod_var,
          " = at_values[[mod_var]], ",
          reg_var,
          " = at_values[[reg_var]])"
        )
      ))
  }
  
  names(mfx_vars) <- NULL

  ## Set up cluster fNULL## Set up cluster for parallelized computation across imputations
  cl <- makeCluster(length(soep_imp$imputations),
                    outfile = "est/stayer_fe_analyses.txt",
                    type = "FORK")
  
  
  
  ## Marginal effects
  mfx <- parLapply(cl,
                   seq_along(est),
                   function (m) {
                     marginaleffects::avg_slopes(
                       model = est[[m]],
                       variables = mfx_vars,
                       by = c(mod_var, reg_var),
                       newdata = pred_data
                     ) %>%
                       as_tibble() %>%
                       dplyr::rename(
                         AME = estimate,
                         SE = std.error,
                         factor = term
                       ) %>%
                       dplyr::select(
                         -contrast,
                         -statistic,
                         -p.value,
                         -s.value,
                         -starts_with("conf."),
                         -starts_with("predicted")
                       )
                   })
  
  ## Exit parallel computation
  stopCluster(cl)
  
  ## Export output
  # Model matrix (stacked)
  mod_summary <- lapply(est, summary)
  
  # Estimates
  data.frame(
    model_id = analysis,
    outcome = model_formula_list_main[[analysis]]$outcome,
    predictor = model_formula_list_main[[analysis]]$predictor,
    moderator = model_formula_list_main[[analysis]]$moderator,
    regional =  model_formula_list_main[[analysis]]$regional,
    owner = model_formula_list_main[[analysis]]$owner,
    spec = model_formula_list_main[[analysis]]$subset,
    rent_var = model_formula_list_main[[analysis]]$rent_var,
    controls = model_formula_list_main[[analysis]]$controls,
    N_obs = nrow(dat[[1]]),
    N_id =  length(unique(dat[[1]]$id)),
    N_plz = length(unique(dat[[1]]$plz)),
    N_kkz = length(unique(dat[[1]]$kkz))
  ) %>%
    write.csv(paste0(
      "est/stayer_fe_analyses/",
      "model_info",
      analysis,
      "_main.csv"
    ))
  
  model_coefficients <- cbind(
    data.frame(
      model_id = analysis,
      x_names = names(coef(est[[1]])),
      b_imp1 = coef(est[[1]]),
      b_imp2 = coef(est[[2]]),
      b_imp3 = coef(est[[3]]),
      b_imp4 = coef(est[[4]]),
      b_imp5 = coef(est[[5]])
    ),
    as.data.frame(as.matrix(vcov(est[[1]]))) %>%
      rename_all(function (x)
        paste0(x, "_imp1")),
    as.data.frame(as.matrix(vcov(est[[2]]))) %>%
      rename_all(function (x)
        paste0(x, "_imp2")),
    as.data.frame(as.matrix(vcov(est[[3]]))) %>%
      rename_all(function (x)
        paste0(x, "_imp3")),
    as.data.frame(as.matrix(vcov(est[[4]]))) %>%
      rename_all(function (x)
        paste0(x, "_imp4")),
    as.data.frame(as.matrix(vcov(est[[5]]))) %>%
      rename_all(function (x)
        paste0(x, "_imp5"))
  ) %>%
    write.csv(paste0(
      "est/stayer_fe_analyses/",
      "model_coefficients",
      analysis,
      "_main.csv"
    ))
  
  pool_mfx_summary(mfx) %>%
    as_tibble() %>%
    mutate(model_id = analysis) %>%
    relocate(model_id, .before = factor) %>%
    write.csv(paste0(
      "est/stayer_fe_analyses/",
      "marginal_effects",
      analysis,
      "_main.csv"
    ))
  
  ## Clear memory
  rm(est, mod_summary, mfx, dat)
  gc(reset = TRUE)
  
  ## Progress
  print(paste0(
    "Analysis: ",
    analysis,
    ". Time elapsed: ",
    round((proc.time()[3] - ptm) / 60, 2),
    " minutes."
  ))
}


## ---- Collect output ----
model_files <- list.files("est/stayer_fe_analyses")
model_info_files_main <-
  model_files[startsWith(model_files, "model_info") &
                endsWith(model_files, "_main.csv")]
model_coefficients_files_main <-
  model_files[startsWith(model_files, "model_coefficients") &
                endsWith(model_files, "_main.csv")]
marginal_effects_files_main <-
  model_files[startsWith(model_files, "marginal_effects") &
                endsWith(model_files, "_main.csv")]
descriptives_files_main <-
  model_files[startsWith(model_files, "descriptives") &
                endsWith(model_files, "_main.csv")]

## ---- Main analyses ----
## Model info
model_info <- data.frame()
for (file in model_info_files_main) {
  model_info <- model_info %>%
    bind_rows(read.csv(paste0("est/stayer_fe_analyses/",
                              file)))
}

model_info <- model_info %>%
  dplyr::select(-X) %>%
  arrange(model_id)

model_info %>%
  write.csv("est/stayer_fe_analyses_model_info_main.csv")


## Model coefficients
model_coefficients <- data.frame()
for (file in model_coefficients_files_main) {
  model_coefficients <- model_coefficients %>%
    bind_rows(read.csv(paste0("est/stayer_fe_analyses/",
                              file)))
}

model_coefficients <- model_coefficients %>%
  dplyr::select(-X) %>%
  arrange(model_id)

model_coefficients %>%
  write.csv("est/stayer_fe_analyses_model_coefficients_main.csv")

## Marginal effects
marginal_effects <- data.frame()
for (file in marginal_effects_files_main) {
  marginal_effects <- marginal_effects %>%
    bind_rows(read.csv(paste0("est/stayer_fe_analyses/",
                              file)))
}

marginal_effects <- marginal_effects %>%
  dplyr::select(-X) %>%
  arrange(model_id)

marginal_effects %>%
  write.csv("est/stayer_fe_analyses_marginal_effects_main.csv")

## Descriptives
descriptives <- data.frame()
for (file in descriptives_files_main) {
  descriptives <- descriptives %>%
    bind_rows(read.csv(paste0("est/stayer_fe_analyses/",
                              file)) %>%
                mutate(model_id = readr::parse_number(file)))
}

descriptives <- descriptives %>%
  dplyr::select(-X) %>%
  arrange(model_id)

descriptives %>%
  write.csv("est/stayer_fe_analyses_descriptives_main.csv")
